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Introduction 



In recent years unstable periodic orbits have been shown to be an effective tool 
in the description of deterministic dynamical systems of low intrinsic dimension Q, 
in diagnosing deterministic chaos in noisy biological systems [^j, and many other 
applications. The theory has been successfully applied to low dimensional ordinary 
differential equations (deterministic chaos) and linear partial differential equations 
(semiclassical quantization). It is an open question whether the theory has anything 
to say about nonlinear partial differential equations (hydrodynamics, field theory). In 
this paper we show that the periodic orbit theory can be used to describe spatially 
extended systems by applying it to the Kuramoto-Sivashinsky equation^, Bj. 

In what follows we shall refer to a periodic solution as a "cycle", and to the 
closure of the union of all periodic solutions as the "invariant set" . Periodic solutions 
are important because they form the skeleton of the invariant setj^, 0], with cycles 
ordered hierarchically; short cycles give good approximations to the invariant set, 
longer cycles refinements. Errors due to neglecting long cycles can be bounded, and 
for nice hyperbolic systems they fall off exponentially or even superexponentially with 
the cutoff cycle length]?], |1. Furthermore, cycles are structurally robust as for smooth 
flows eigenvalues of short cycles vary slowly with smooth parameter changes, short 
cycles can be accurately extracted from experimental or numerical data, and global 
averages (such as correlation exponents, escape rates and other "thermodynamic" 
averages) can be efficiently computed from short cycles by means of cycle expansions. 

While the role of periodic solutions in elucidating the asymptotics of ordinary 
differential equations was already appreciated by Poincare||, allegedly HopfJTot 
and, more demonstrably, Spiegel and collaborators [jllj [l^, [l3| have argued that the 
asymptotics of partial differential equations should also be described in terms of 
recurrent spatiotemporal patterns. Pictorially, dynamics drives a given spatially 
extended system through a repertoire of unstable patterns; as we watch a given 
"turbulent" system evolve, every so often we catch a glimpse of a familiar pattern. 
For any finite spatial resolution, the system follows approximately for a finite time 
a pattern belonging to a finite alphabet of admissible patterns, and the long term 
dynamics can be thought of as a walk through the space of such patterns, just as 
chaotic dynamics with a low dimensional attractor can be thought of as a succession 
of nearly periodic (but unstable) motions. 



1. Kuramoto-Sivashinsky system 

We offer here a modest implementation of the above program on a prototype spatially 
extended dynamical system defined by the Kuramoto-Sivashinsky equation || [f| 

V>t (y> )x ^xx fUxxxx (1) 

which arises as a model amplitude equation for interfacial instabilities in a variety of 
contexts - see e.g. reference Here t > is the time, x 6 [0, 2ir) is the space 

coordinate, and v is a fourth-order "viscosity" damping parameter. The subscripts x 
and t denote the partial derivatives with respect to x and t. We take the Kuramoto- 
Sivashinsky system because it is one of the simplest physically interesting spatially 
extended nonlinear systems, but in the present context the interpretation of the 
equation, or the equation itself is not the most important ingredient; the approach 
should be applicable to a wide class of spatially extended nonlinear systems. The 



3 



salient feature of such partial differential equations is that for any finite value of 
v their asymptotics is in principle describable by a finite set of "inertial manifold" 
ordinary differential equations |l5|]. 

The program of studying unstable solutions in this context was initiated by Goren, 
Eckmann and Procaccia[[l6] who have used a 2-unstable modes truncation of the 
Kuramoto-Sivashinsky equation to study the dynamics connecting coexisting unstable 
temporally stationary solutions. We shall study here unstable spatiotemporally periodic 
solutions of the full Kuramoto-Sivashinsky system. Our main result is that in the limit 
of weak turbulence or "spatiotemporal chaos" , we can determine hierarchically and 
exhaustively cycles of longer and longer periods, and apply this data to the evaluation 
of global averages. 

The function u(x, f) = u(x + 2-7T, t) is assumed periodic on the x £ [0, 2n] interval. 
As u(x,t) has compact support, the standard strategy is to expand it in a discrete 
spatial Fourier series: 

+00 

u(x,t)= h(t)e ikx . (2) 

k— — oo 

Since u(x,t) is real, bk = b*_ k . Substituting (Q) into (jl|) yields the infinite ladder of 
evolution equations for the Fourier coefficients bk' 

00 

h = (k 2 - vk A )b k + ifc ^ b m b k - m • (3) 

m — — oo 

As bo = 0, the average (the mean drift) of the solution is an integral of motion. In 
what follows we shall assume that this average is zero, J dx u(x, t) = 0. 

The coefficients bk are in general complex functions of time. We can simplify the 
system (|3|) further by assuming that bk are pure imaginary, bk — ia^, where are real. 
As we shall see below, this picks out the subspace of odd solutions u(x, t) — —u(—x, i), 
with the evolution equations 

00 

hk = (k 2 - vk 4 )a k - fc a m ak- m ■ (4) 

m— — 00 

We shall determine the periodic solutions in the space of Fourier coefficients, and then 
reconstitute from them the unstable spatiotemporally periodic solutions of ([!]). 

The trivial solution u(x,t) = is a fixed point of (|l|). From (Q) it follows that 
the |fc| < 1/ ' sjv long wavelength modes of this fixed point are linearly unstable, and 
the |fc| > short wavelength modes are stable. For v > 1, u(x.t) — is the 

globally attractive stable fixed point; starting with v = 1 the solutions go through a 



rich sequence of bifurcations, studied e.g. in reference |14 . Detailed knowledge of the 
parameter dependence of bifurcations sequences is not needed for our purposes; we 
shall take sufficiently small so that the dynamics can be spatiotemporally chaotic, 
but not so small that we would be overwhelmed by too many short wavelength modes 
needed in order to accurately represent the dynamics. 

The growth of the unstable long wavelengths (low |fc|) excites the short wavelengths 
through the nonlinear term in (^). The excitations thus transferred are dissipated 
by the strongly damped short wavelengths, and a sort of "chaotic equilibrium" can 
emerge. The very short wavelengths |fc| 3> wm remain small for all times, 

but the intermediate wavelengths of order ~ wm pl a Y an important role 
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in maintaining the dynamical equilibrium. As the clamping parameter decreases, the 
solutions increasingly take on Burgers type shock front character which is poorly 
represented by the Fourier basis, and many higher harmonics need to be kept Jl4|, |l6| 
in truncations of (Q). Hence, while one may truncate the high modes in the expansion 
, care has to be exercised to ensure that no modes essential to the dynamics are 
chopped away. 

Before proceeding with the calculations, we take into account the symmetries of 
the solutions and describe our criterion for reliable truncations of the infinite ladder 
of ordinary differential equations (0). 



2. Symmetry decomposition 

As usual, the first step in analysis of such dynamical flows is to restrict the dynamics to 
a Poincare section. We shall fix the Poincare section to be the hyperplane a\ — 0. We 
integrate (Q) with the initial conditions a\ = 0, and arbitrary values of the coordinates 
ci2, . . . , aw, where N is the truncation order. When a\ becomes the next time, the 
coordinates a^, ■ ■ ■ , % are mapped into (a' 2 , . . . a' N ) — P(a 2 , ■ ■ ■ , ctj\r), where P is the 
Poincare map. P defines a mapping of a N — 1 dimensional hyperplane into itself. 
Under successive iterations of P, any trajectory approaches the attractor A, which 
itself is an invariant set under P. 

A trajectory of (^) can cross the plane a\ — in two possible ways: either when 
a'i > ("up" intersection) or when a'i < ("down" intersection), with the "down" 
and "up" crossings alternating. It then makes sense to define the Poincare map P as 
a transition between, say, "up" and "up" crossing. With Poincare section defined as 
the "up-up" transition, it is natural to define a "down-up" transition map 0. Since 
describes the transition from down to up (or up to down) state, the map 2 describes 
the transition up-down-up, that is 2 = P. 

Consider the spatial flip and shift symmetry operations Ru(x) = u(—x), Su(x) — 
u(x + 7r). The latter symmetry reflects the invariance under the shift u(x, t) — > 
u(x + 7T, t), and is a particular case of the translational invariance of the Kuramoto- 
Sivashinsky equation (Q) . In the Fourier modes decomposition (^) this symmetry acts 
as S : d2k — > <22fc)fl2/c+i - * —a-ik+i- Relations R 2 = S 2 = 1 induce decomposition of 
the space of solutions into 4 invariant subspacesjOJ; the above restriction to bk = iflfc 
amounts to specializing to a subspace of odd solutions u(x, t) = —u(—x, t). 

Now, with the help of the symmetry S the whole attractor Atot can be decomposed 
into two pieces: Atot = Ao U SAo for some set _4o- It can happen that the set Ao 
(the symmetrically decomposed attractor) can be decomposed even further by the 
action of the map 0. In this case the attractor will consist of four disjoint sets: 
Atot = «4U5^4U0^4U05^4. As we shall see, this decomposition is not always possible, 
since sometimes A overlaps with 05^1 (in this case 0.4 will also overlap with SA). 
We shall carry out our calculations in the regime where the decomposition into four 
disjoint pieces is possible. In this case the set A can be taken as the fundamental 
domain of the Poincare map, with 5^4, 0.4 and 05^1 its images under the S and 
mappings. 

This reduction of the dynamics to the fundamental domain is particularly useful 
in periodic orbit calculations, as it simplifies symbolic dynamics and improves the 
convergence of cycle expansions ]l7[. 
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3. Fourier modes truncations 

When we simulate the equation (||) on a computer, we have to truncate the ladder of 
equations to a finite length N, i.e., set cik — for fc > N. N has to be sufficiently 
large that no harmonics important for the dynamics with k > N are truncated. 
On the other hand, computation time increases dramatically with the increase of iV: 
since we will be evaluating the stability matrices for the flow, the computation time 
will grow at least as iV 2 . 

Adding an extra dimension to a truncation of the system ([|) introduces a small 
perturbation, and this can (and often will) throw the system into a totally different 
asymptotic state. A chaotic attractor for N = 15 can become a period three window 
for N = 16, and so on. If we compute, for example, the Lyapunov exponent A(V, N) 
for the strange attractor of the system (||), there is no reason to expect \(v,N) to 
smoothly converge to the limit value X(v,oo) as N — > oo. The situation is different 
in the periodic windows, where the system is structurally stable, and it makes sense 
to compute Lyapunov exponents, escape rates, etc. for the repeller, i.e. the closure of 
the set of all unstable periodic orbits. Here the power of cycle expansions comes in: 
to compute quantities on the repeller by direct averaging methods is generally more 
difficult, because the motion quickly collapses to the stable cycle. 

We have found that the minimum value of N to get any chaotic behavior at all was 
N = 9. However, the dynamics for the N = 9 truncated system is rather different from 
the full system dynamics, and therefore we have performed our numerical calculations 
for N = 15, N = 16 and N = 17. Figure [l] is a representative plot of the Feigenbaum 
tree for the Poincare map P. To obtain this figure, we took a random initial point, 
iterated it for a some time to let it settle on the attractor and then plotted the 
coordinate of the next 1000 intersections with the Poincare section. Repeating this for 
different values of the damping parameter v, one can obtain a picture of the attractor 
as a function of v. For an intermediate range of values of v, the dynamics exhibits a 
rich variety of behaviours, such as strange attractors, stable limit cycles, and so on. 
The Feigenbaum trees for different values of N resemble each other, but the precise 
values of v corresponding the various bifurcations depend on the order of truncation 
N. 

Based on the observed numerical similarity between the Feigenbaum trees for 
N = 16 and N — 17 (cf. figure Q), we choose N — 16 as a reasonable cutoff and 
will use only this truncation throughout the remainder of this paper. We will examine 
two values of the damping parameter: v = 0.029910, for which the system is chaotic, 
and v — 0.029924, for which the system has a stable period-3 cycle. In our numerical 
work we use both the pseudospectral|ll| and the 4-th order variable-step Runge-Kutta 
integration routines |Q; their results are in satisfactory agreement. As will be seen 
below, the good control of symbolic dynamics guarantees that we do not miss any short 
periodic orbits generated by the bifurcation sequence indicated by the Feigenbaum tree 
of figure [|. However, even though we are fairly sure that for this parameter value we 
have all short periodic orbits, the possibility that other sets of periodic solutions exist 
somewhere else in the phase space has not been excluded. 

The problem with such high dimensional truncations of (|J) is that the dynamics 
is difficult to visualize. We can examine its projections onto any three axes 
ai,a,j,a,k, as in figure or, alternatively, study a return map for a given coordinate 
afc — > a' k = Pk{d2, ■ ■ ■ ,o-n) as the one plotted in figure The full return map is 
[N — l)-dimensional a — > P(a2, . . . ,ajv) = a' and single-valued, and for the values 
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Sgure 1. Fcigcnbaum tree for coordinate a§, TV = 16 Fourier modes truncation of 
). The two upper arrows indicate the values of damping parameter that we use 
in our numerical investigations; v = 0.029910 (chaotic) and v = 0.029924 (period-3 
window). The lower arrow indicates the kink where the invariant set A starts to 
overlap with OSA. Truncation to TV = 17 modes yields a similar figure, with values 
for specific bifurcation points shifted by ~ 10 -5 with respect to the TV = 16 values. 
The choice of the coordinate aa is arbitrary; projected down to any coordinate, the 
tree is qualitatively the same. 

of v used here the attractor is essentially 1-dimensional, but its projection into the 
{at, Pk(a>2, ■ ■ ■ , ajv)} plane can be multi- valued and self-intersecting. One can imagine 
a situation where no "good" projection is possible, that is, any projection onto any 
two-dimensional plane is a multiple- valued function. The question is how to treat such 
a map? 

4. One-dimensional visualization of the dynamics 

We now describe an approach which simplifies matters a lot by reducing the map to an 
approximate one-dimensional map. The multiple- valuedness in figure ^ arises from the 
fact that the return map is a 2-dimensional projection of a convoluted 1-dimensional 
curve embedded into a high-dimensional space. We shall show that it is possible to 
find an intrinsic parametrization s along the unstable manifold, such that the map 
s — > f(s) induced by the full d-dimensional flow is approximately 1-dimensional. 
Strictly speaking, the attractor on figure has a certain thickness transverse to it, 
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Figure 2. Projections of a typical 16-dimensional trajectory onto different 3- 
dimensional subspaces, coordinates (a) {01,02,03}, (b) {01,02,04}. N = 16 Fourier 
modes truncation with v = 0.029910. 
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Figure 3. The attractor of the system (^), plotted as the ae component of the 01 = 
Poincare section return map, 10,000 Poincare section returns of a typical trajectory. 
Indicated are the periodic points 0, 1 and 01; as this is an arbitrary projection of 
the invariant set, they exhibit no good spatial ordering. TV = 16 Fourier modes 
truncation with v = 0.029910. 
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but the contraction in the transverse directions is so strong that the invariant set is 
effectively 1-dimensional. 

Suppose we already have determined some of the shorter cycles for our system, i.e. 
the fixed points of the Poincare map and its iterates. This is accomplished relatively 
easily by checking a trajectory of a random initial point for close returns and then 
using these as initial guesses for a cycle search algorithm. We now assume that the 
invariant set can be approximated by a curve passing close to all periodic points, and 
determine the order of periodic points along such curve. This is done in the following 
way: there exists a fixed point which is not connected to the attractor (the point 
in figure ||) - we choose this fixed point as the starting point and assign it number 
1. Point number 2 is the periodic point in the sample which is closest (in the full 
space) to this fixed point, and the n-th point is determined as the point which has the 
minimum distance from the point number n — 1 among all the periodic points which 
have not yet been enumerated. Proceeding this way, we order all the periodic points 
that we have found so far. 

Since all periodic points belong to cycles, their images are known and are simply 
the successive periodic points along the cycle. We use this fact to recursively construct 
a 1-dimensional mapping Si — > f(si). We approximate parametrization length s along 
the invariant set by computing the Euclidean inter-distances between the successive 
periodic points in the full dynamical space, s\ = 0, S2 = \\bl2 — a_i || , — s,;_i = 
||aj — SLj_i || . The i-th cycle point Si is mapped onto its image s a i — f(si), where ai 
denotes the label of the next periodic point in the cycle. We can now find longer 
periodic orbits of the 1-dimensional map / by standard methods such as inverse 
iteration, and guess the location of the corresponding points in the full iV-dimensional 
space by interpolating between the nearest known periodic points. These will not be 
exact periodic orbits of the full system, but are very useful as good starting guesses in 
a search for the exact periodic orbits. Iteratively, more and more periodic orbits can 
be computed. While it only pays to refine the 1-dimcnsional parametrization until the 
density of the periodic points become so high that the width of the attractor becomes 
noticeable, the 1-dimensional map continues to provide good initial guesses to longer 
periodic orbits. More sophisticated methods are needed only if high accuracy around 
the folding region of /(s) is required in order to distinguish between long cycles. 

For the values of v we are working with, the attractor consists of four disjoint 
sets, the fundamental domain A and its images under the maps S and 0. In this case 
the approximate return map s — > f(s) is unimodal. The corresponding map on the 
symmetric part of the attractor, SOA, is likewise unimodal, and turned 180 degrees 
around the origin. For the values of v we work with the two maps do not interact 
and their domains are separate. However, if the value of the damping parameter v is 
decreased sufficiently, the domains of the maps join and together they form a connected 
invariant set described by a bimodal map p0| . This joining of the fundamental domain 
A and its symmetric image OS A is visible in figure [l] at v ~ 0.0299, where the range 
of the clq coordinate increases discontinuously. 

We use the unimodal map s — ► /(s) to construct binary symbolic dynamics for the 
system in the usual way: assign the symbol '0' to points to the left of the maximum, 
T' to the points to the right. In the period-3 window with the stable cycle 001, the 
pruning rules are very easy: except for the stable 001 cycle and the fixed point 
(both disjoint from the invariant set) two 0's in a row are forbidden. In this case 
it is convenient to redefine the alphabet by denoting the symbol pair 01 by a and 
the symbol 1 by b. This renders the symbolic dynamics of the points on the repeller 
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Figure 4. The return map s n _|_i = f(s n ) constructed from the images of periodic 
points. The diamonds were obtained by using 34 periodic points, the dots were 
obtained by using 240 periodic points. We have indicated the periodic points 0, 1 
and 01. Note that the transverse fractal structure of the map shows when the number 
of points is increased. N = 16 Fourier modes truncation with v = 0.029910. 
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Figure 5. The return timeJT(s) as a function of the parameter s, evaluated on the 
periodic points, as in figure W with the diamonds obtained by 34 periodic points and 
the dots by 240 periodic points. The fine structure is due to the fractal structure of 
the attractor. 



10 



complete binary: all sequences of the letters a and b are admissible. 

A flow in N dimensions can be reduced to an (N — l)-dimensional return map 
by suspension on a Poincare section provided that the Poincare return map is 
supplemented by a "time ceiling function" |2lj] which accounts for a variation in the 
section return times. Hence we also determine the return time T(si) for each periodic 
point i, and use those to construct recursively the periodic orbit approximations to 
the time ceiling function, figure §. The mean Poincare section return time is of order 
T « .88. 

4-.1. Numerical results 

We have found all cycles up to topological length 10 (the redefined topological length in 
the case of the period-3 window), 92 cycles in the chaotic regime and 228 in the period- 
3 window, by using the 1-dimensional parametrization /(s) to find initial guesses for 
periodic points of the full N — 16 Fourier modes truncation and then determining the 
cycles by a multi-shooting Newton routine. It is worth noting that the effectiveness 
of using the 1-dimensional f(s) approximation to the dynamics to determine initial 
guesses is such that for a typical cycle it takes only 2-3 Newton iterations to find the 
cycle with an accuracy of 10~ 10 . 
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Figure 6. Lyapunov exponents Aj. versus k for the periodic orbit 1 compared with 
the stability eigenvalues of the u(x,t) = stationary solution k 2 — uk 4 . \k for k > 8 
fall below the numerical accuracy of integration and are not meaningful. N = 16 
Fourier modes, v = 0.029924, chaotic regime. 

In table |l| we list the periodic orbits to topological length 5 found by our method. 
The value of A2 serves as an indication of the accuracy of our numerics, as A2 
corresponds to the marginal eigenvalue along the periodic orbit, strictly equal to 1. 
All cycles seem to have real eigenvalues (to within the numerical accuracy) except 
for the 0-cycle which has a pair of complex eigenvalues, A3 and A4. We therefore 
do not list the corresponding imaginary parts of the eigenvalues. To illustrate the 
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Table 1. All cycles up to topological length 5 for the N = 16 Fourier modes 
truncation of the Kuramoto-Sivashinsky equation (^), damping parameter v = 
0.029910 (chaotic attractor) and v = 0.029924 (period-3 window), their itineraries, 
periods, the first four stability eigenvalues. For the chaotic attractor pruning shows 
up at the topological length 4; 0001 and 0011 cycles are pruned. The deviation from 
unity of A2, the eigenvalue along the flow, is an indication of the accuracy of the 
numerical integration. For the period-3 window we also give the itineraries in the 
redefined alphabet where a = 1 and b = 10. 



p 




T p 


Ai 


A 


2 — 


1 


A3 






A 4 






Chaotic, v = 


0.029910 



























0.897653 


3.298183 


5 


10" 


12 


-2.793085 


10" 


3 


-2.793085 


10" 


3 


1 




0.870729 


-2.014326 


-5 


10~ 


12 


6.579608 


10~ 


3 


-3.653655 


10" 


4 


10 




1.751810 


-3.801854 


8 


10~ 


12 


-3.892045 


10" 


5 


2.576621 


10- 


7 


100 




2.639954 


-4.852486 


1 


10" 


11 


3.044730 


10" 


7 


-3.297996 


10- 


10 


110 




2.632544 


6.062332 


2 


10" 


11 


-2.721273 


10" 


7 


-1.961928 


10- 


10 


1000 


























1100 


























1110 




3.497622 


-14.76756 


2 


10" 


11 


-1.629532 


10" 


9 


6.041192 


10" 


14 


10100 




4.393973 


19.64397 


2 


10" 


11 


-1.083266 


10" 


11 


3.796396 


10- 


15 


11100 




4.391976 


-18.93979 


2 


10" 


11 


1.162713 


10" 


11 


-1.247149 


10" 


14 


11010 




4.380100 


-26.11626 


2 


10" 


11 


1.005397 


10" 


11 


8.161650 


10- 


15 


11110 




4.370895 


28.53133 


2 


10" 


11 


1.706568 


10" 


11 


1.706568 


10" 


14 


Period-3 window, v = 0.029924 







0.897809 


3.185997 


7 


10" 


13 


-2.772435 


10" 


3 


-2.772435 


10- 


3 


1 


a 


0.871737 


-1.914257 


5 


10" 


13 


6.913449 


10" 


3 


-3.676167 


10" 


4 


10 


b 


1.752821 


-3.250080 


1 


10" 


12 


-4.563478 


10~ 


5 


2.468647 


10" 


7 


100 




2.638794 


-0.315134 


-4 


10" 


13 


4.821809 


10" 


6 


-2.576341 


10" 


10 


110 


ab 


2.636903 


2.263744 


3 


10" 


12 


-6.923648 


10- 


7 


-2.251226 


10- 


10 


1110 


aab 


3.500743 


-10.87103 


2 


10" 


12 


-2.198314 


10- 


9 


3.302367 


10- 


14 


11010 


abb 


4.382927 


-15.84102 


2 


10" 


12 


1.656690 


10- 


11 


1.388232 


10- 


14 


11110 


aaab 


4.375712 


18.52766 


3 


10" 


12 


-1.604898 


10" 


11 


2.831886 


10" 


14 



rapid contraction in the nonleading eigendirections we plot all the eigenvalues of the 
1-cycle in figure ^. As the length of the orbit increases, the magnitude of contracting 
eigenvalues falls quickly bellow the attainable numerical numerical accuracy « 10 -16 
and our numerical results for A& are not meaningful for k > 8. 

Having determined the periodic solutions p in the Fourier modes space, we now 
go back to the configuration space and plot the corresponding spatiotemporally 
periodic solutions u p (x,t): they are the repertoire of the recurrent spatiotemporal 
patterns that Hopf wanted to see in turbulent dynamics. Different spatiotemporally 
periodic solutions are qualitatively very much alike but still different, as a closer 
inspection reveals. In figure [?] we plot uo(x,t) corresponding to the Fourier space 
0-cycle. Other solutions, plotted in the configuration space, exhibit the same overall 
gross structure. For this reason we find it more informative to plot the difference 
uo(x,t'To) — u p {x,t"T p /n p ) rather than u p {x,t) itself. Here p labels a given prime 
(non-repeating) cycle, n p is the topological cycle length, T p its period, and the time 
is rescaled to make this difference periodic in time: t' — t/To and t" — n p t/T p , 
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so that t" ranges from to n p . u (x,t'To) — Ui(x,t"T 1 ) is given in figure |], and 
u Q (x,t'T ) - u i(x,t"T 01 /2) in figure [| 




Figure 7. Spatiotcmporally periodic solution uo(x,t). We have divided x by it and 
plotted only the x > part, since we work in the subspace of the odd solutions, 
u(x,t) = —u(—x,t). N = 16 Fourier modes truncation with u = 0.029910. 




Figure 8. The difference between the two shortest period spatiotcmporally periodic 
solutions uo(x,t'To) and u\(x,t"T\). 



5. Global averaging: periodic orbits in action 

The above investigation of the Kuramoto-Sivashinsky system demonstrates that it 
is possible to construct recursively and exhaustively a hierarchy of spatiotemporally 
periodic unstable solutions of a spatially extended nonlinear system. 

Now we turn to the central issue of this paper; qualitatively, these solutions are 
indeed an implementation of Hopf's program, but how is this information to be 
used quantitatively? This is precisely what the periodic orbit theory is about; it 
offers machinery that puts together the topological and the quantitative information 
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Figure 9. The difference between solution uo(x,t'To) repeated twice and the n p = 2 
period spatiotemporally periodic solution uoi(x,t"Toi/2). 



about individual solutions, such as their periods and stabilities, into predictions about 
measurable global averages, such as the Lyapunov exponents, correlation functions, 
and so on. The proper tool for computing such global characterizations of the dynamics 
are the trace and determinant formulas of the periodic orbit theory. 

We shall briefly summarize the aspects of the periodic orbit theory relevant to the 
present application; for a complete exposition of the theory the reader is referred to 
reference [^2|. The key idea is to replace a time average &(x)/t of an "observable" <fi 
measured along a dynamical trajectory x(t) — f l (x) 

&( X ) = [ dT(f>(x(T)) 



by the spatial average ^e' 3 * J of the quantity e' 3 * ^ . Here (3 is a dummy variable, used 
to recover the desired expectation value {4>) = lim^oo (3>*/t) by taking derivatives 
of (e' 3 ** ) and then setting j3 = 0. For large t the average ( e' 3 *' ) behaves as the trace 



00 e rf3<S> p 

^'°£r p £ l ^-rr,,) (5) 

p r—1 I V VI I 

of the evolution operator 

£ t (x,y) = 6(y-f t (x))e^ t ^. 

and is dominated by its largest eigenvalue e ts ^\ 

The trace formula (|^) has an intuitive geometrical interpretation. The sums in 
(|J) are over prime periodic orbits p and their repeats r, T p are their periods, and Jp 
are their stability matrices. Prime cycles partition the dynamical space into closed 
tubes of length T p and thickness |det(l — J p )| , and the trace picks up a periodic 
orbit contribution only when the time t equals a prime period or its repeat, hence the 
time delta function 5(t — rT p ). Finally, e /3 * p is the mean value of e^ $ ^ evaluated on 

this part of dynamical space, so the trace formula is the average of (e' 3 * ) expressed 
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as a partition of the space of solutions into a repertoire of spatiotemporally periodic 
solutions, each weighted by its stability, i.e. likelihood of its occurrence in a long time 
evolution of the system. 

In applications of the periodic orbit theory the related Fredholm determinant 




has better convergence as a function of the maximal cycle length truncation, so that is 
the function whose leading zero F((3, s{f3)) = we determine here in order to evaluate 
the leading eigenvalue s(/3). 

The dummy variable z in (^|) keeps track of the topological lengths n p (number of 
the Poincare section crossings), and is used to expand F as a series in z. If we know 
all cycles up to topological length I we truncate F to Z-th order polynomial: 

F({3, s) = 1 - c k z k + (remainder) (7) 
l 

and set z = 1. The general theory [^3[ ^[ |) then guarantees that for a hyperbolic 
dynamical system the coefficients fall off in magnitude exponentially or faster with 
increasing k. We now calculate the leading eigenvalue s(/3) by determining the smallest 
zero of F{(3, s), and check the convergence of this estimate by studying it as a function 
of the maximal cycle length truncation I. If the flow conserves all trajectories, the 
leading eigenvalue must satisfy s(0) = 0; if the invariant set is repelling, the leading 
eigenvalue yields 7 = — s(0), the escape rate from the repeller. Once the leading 
eigenvalue is determined we can calculate the desired average (</)) using formulaH: 



<*> - -I 



dF JdF 



(8) 

»=»(o) 



For example, if we take as our "observable" log |A* (x)|, the largest eigenvalue of the 
linearized stability of the flow, <E> P will be log |Ai p | where Ai )P is the largest eigenvalue 
of stability matrix of the cycle p, and the above formula yields the Lyapunov exponent 
(A). 

Both the numerator and the denominator in (|J) have a cycle expansion analogous 
to (0) (cf. refer ence j22j), and the same periodic orbit data suffices for their evaluation. 

Conceptually the most important lesson of the periodic orbit theory is that the 
spatiotemporally periodic solutions are not to be thought of as eigenmodes to be used 
as a linear basis for expressing solutions of the equations of motion - as the equations 
are nonlinear, the periodic solutions are in no sense additive. Nevertheless, the trace 
formulas and determinants of the periodic orbit theory give a precise prescription for 
how to systematically explore the repertoire of admissible spatiotemporal patterns, 
and how to put them together in order to predict measurable observables. 



5.1, Numerical results 

One of the objectives of a theory of turbulence is to predict measurable global averages 
over turbulent flows, such as velocity-velocity correlations and transport coefficients. 
While in principle the periodic orbit averaging formulas should be applicable to such 
averages, with the present parameter values we are far from any strongly turbulent 
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Figure 10. log lc) of the coefficients |cfe| in the cycle expansion (?| of F(0, 0) versus 
k for the period-3 window case (crosses) and the chaotic case (diamonds). N = 16 
Fourier modes truncation. 



regime, and here we shall restrict ourselves to the simplest tests of chaotic dynamics: 
we shall test the theory by evaluating Lyapunov exponents and escape rates. 

We compute the periodic orbits, escape rates and Lyapunov exponents both for the 
period-3 window and a chaotic regime. In the case of period-3 window the complete 
symbolics dynamics and grammar rules are known and good convergence of cycle 
expansions is expected both for the escape rate from the repeller and the Lyapunov 
exponent. Parenthetically, the stable period-3 orbit is separated from the rest of the 
invariant set by its immediate basin of attraction window, and its eigenvalues bear no 
immediate relation to the escape rate and the Lyapunov exponent of the repelling set. 

In the case of a generic "strange attractor", the convergence is not expected to 
be nearly as good, since in this case there exist no finite description of the symbolic 
dynamics. For closed systems (no escape) 7 = and F (0, 0) = 0. The discrepancy 
of the value ^(0, 0) from for a closed system allows us to estimate the accuracy of 
finite cycle length approximations to the Fredholm determinant. 

The analytic properties of the Fredholm determinant are illustrated by the decay 
rate of the coefficients Ck as a function of k in the expansion ((7J). If the complete 
symbolic dynamics is known and the system is hyperbolic, the decay of Ck should be 
superexponential[Q. This is illustrated in figure [To], where we plot the coefficients Ck 
for the 16-dimensional system for the chaotic case and for the period-3 window. We 
can clearly see the superexponential decay for the period-3 window case and at best 
exponential decay for the chaotic case. 

Our results are presented in table ||. One observes that when the symbolic 
dynamics is known (period-3 window), the convergence is much better than in the 
generic case, in accordance with the periodic orbit theory expectations. 
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Table 2 . The escape rate 7 and the leading Lyapunov exponent as a function of the 
cycle expansion truncation n max for the N = 16 Fourier modes truncation, chaotic 
regime (y = 0.029910) and period-3 window [y = 0.029924). In the period-3 window 
the Fredholm determinant starts converging only for n m ax > 4; for Umax = 4 it has 
no real zero at all. A numerical simulation estimate for the Lyapunov exponent in 
the chaotic regime is given in the last line; for this parameter value the escape rate, 



'Yj should 


strictly ecjual 










chaotic 




period-3 


window 




7 


Ai 


7 


Ai 


1 






0.428143 


0.703010 


2 


0.441948 


0.981267 


-0.187882 


0.430485 


3 


0.080117 


0.765050 


-0.049325 


0.469350 


4 


0.148583 


0.703072 






5 


0.068513 


0.727498 


1.072468 


0.585506 


6 


0.027724 


0.699907 


0.078008 


0.547005 


7 


0.035137 


0.693852 


0.088132 


0.598977 


8 


0.007104 


0.675529 


0.090425 


0.631551 


9 


0.021066 


0.673144 


0.090101 


0.618160 


10 


0.007367 


0.646233 


0.090065 


0.621271 


numer. 




0.629 







6. Summary 

Hopf's proposal for a theory of turbulence was, as we understand it, to think of 
turbulence as a sequence of near recurrences of a repertoire of unstable spatiotemporal 
patterns. Hopf's proposal is in its spirit very different from most ideas that animate 
current turbulence research. It is distinct from the Landau quasiperiodic picture 
of turbulence as a sum of infinite number of incommensurate frequencies, with 
dynamics taking place on a large-dimensional torus. It is not the Kolmogorov's 
1941 homogeneous turbulence with no coherent structures fixing the length scale, 
here all the action is in specific coherent structures. It is emphatically not universal; 
spatiotemporally periodic solutions are specific to the particular set of equations 
and boundary conditions. And it is not probabilistic; everything is fixed by the 
deterministic dynamics with no probabilistic assumptions on the velocity distributions 
or external stochastic forcing. 

Our investigation of the Kuramoto-Sivashinsky system is a step in the direction 
of implementing Hopf's program. We have constructed a complete and exhaustive 
hierarchy of spatiotemporally periodic solutions of spatially extended nonlinear system 
and applied the periodic orbit theory to evaluation of global averages for such 
system. Conceptually the most important lesson of this theory is that the unstable 
spatiotemporally periodic solutions serve to explore systematically the repertoire of 
admissible spatiotemporal patterns, with the trace and determinant formulas and their 
cycle expansions being the proper tools for extraction of quantitative predictions from 
the periodic orbits data. 

We have applied the theory to a low dimensional attractor, not larger than the 
Lorenz's original strange attractor [p4|. As our aim was to solve the given equations 
accurately, we were forced to work with a high dimensional Fourier modes truncations, 
and we succeeded in determining the periodic orbits for flows of much higher dimension 
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than in previous applications of the periodic orbit theory. As something new, we have 
developed an intrinsic parametrization of the invariant set that provided the key to 
finding the periodic orbits. 

In practice, the method of averaging by means of periodic orbits produced best 
results when the complete symbolic dynamics was known. For generic parameter 
values we cannot claim that the periodic orbit approach is computationally superior 
to a direct numerical simulation. A program to find periodic orbits up to length 10 for 
one value of the damping parameter v requires a day of CPU on a fast workstation, 
much longer than the time used in the direct numerical simulations. 

The parameter v values that we work with correspond to the weakest nontrivial 
"turbulence" , and it is an open question to what extent the approach remains 
implementable as the system goes more turbulent. Our hope is that the unstable 
structures captured so far can be adiabatically tracked to the "intermediate 
turbulence" regime, and still remain sufficiently representative of the space of 
admissible patterns to allow meaningful estimates of global averages. As long as no 
effective coordinatization of the "inertial manifold" exists and we rely on the spatial 
Fourier decomposition, the present approach is bound to fail in the "strong turbulence" 
v — ► limit, where the dominant structures are Burgers- type shocks and truncations 
of the spatial Fourier modes expansions are increasingly uncontrollable. 
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